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Abstract 

This paper discusses aspects of the second order hyperboUc partial differential equation associ- 
ated with the ideal lossless string under tension and it's relationship to two discrete models. These 
models are finite differencing in the time domain and digital waveguide models. It is known from 
the theory of partial differential operators that in general one has to expect the string to accumulate 
displacement as response to impulsive excitations. Discrete models should be expected to display 
comparable behavior. As a result it is shown that impulsive propagations can be interpreted as 
the difference of step functions and hence how the impulsive response can be seen as one case of 
the general integrating behavior of the string. Impulsive propagations come about in situations of 
time-symmetry whereas step-function occur as a result of time-asymmetry. The difference between 
the physical stability of the wave equation, which allows for unbounded growth in displacement, 
and computational stability, that requires bounded growth, is derived. 
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I. INTRODUCTION 



Our purpose here is to discuss aspects of the relationship of the solution of the one- 
dimensional second order wave equation to two discrete models thereof. The first discrete 
model is the digital waveguide model in one spatial dimension. The second discrete model 
is a finite difference model in the time domain. In particular we will also discuss how this 



relationship explains different behavior 



between the discrete models. This relationship has 



drawn much attention recentlyjl|, O, S, I4 , 0, 0, Q, 0, 1^ 

It is shown that the finite difference model can account for solutions of the wave equation 
and that these solutions are physically meaningful. 

This allows for a direct interpretation of recent results by Karjalainen and Erkut ^ from 
the fundamental solution of the wave equation. Karjalainen and Erkut gave the restrict- 
ing conditions necessary to make finite difference models in the time domain and digital 
waveguide models connect. 

Regarding the stability behavior of the discrete models, the continuous stability is dis- 
cussed and it is shown that physically stable responses of the string may appear unstable in 
a discrete model or signal-processing sense. 

The paper is structured as follows. First known derivations of the solution of the wave 

utions of its partial 



1^ 



and the digital 



equation is given, both classically and via the theory of fundamental so 
differential operators. A discrete comparison of the finite difference 
waveguide model j3] follows. Smith's textj^ provides the authoritative summary in of digital 
waveguides with respect to the continuous derivations given earlier. We conclude with 
implications of these observations. 



II. SOLUTION OF THE WAVE EQUATION IN ONE DIMENSIONS 

The results in this section are well-known. They are repeated here to facilitate arguments 
in the following sections. Of concern is the general digital simulation of a string under force. 
The free ideal string is well-described by the the 1 + 1 dimensional d'Alembertian operator 
on a scalar field I 



OyM = i^2-^'§2)y(^^t) (1) 
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An external force leads to the inhomogeneous case of the same equation: 



ny{x,t) = fix,t) (2) 

Without loss of generality we assume that c = 1 for our discussion and hence one gets 
the factored form of the d'Alembertian: 



~ \di^ Wt) \di~ di) ^' 

The factored form suggests the substitution ^ = x — t and rj = x + t we obtain the 
canonical form of the wave equation: 



d^dr] 4' V 2 ' 2 



This can be directly integrated, yielding 







rv /"€ 

u= / (j){u,w) dudw + hi{^) + h2{r]) (5) 
•J no -^io 

where hi, and h2 are "constants of integration". However, this notation hides that either 
one of these constants has been integrated over with regards to its parameter. Also these 
are not uniquely defined functions from a yet undefined functional space but rather any 
function from suitable family of functions The derivation procedure suggests that /ii(-) 
and /i2(') are one and twice differentiable everywhere in the solution space, but which one 
is twice differentiable depends in what order the solution has been integrated over. 

The solution of the homogeneous case (^J) with initial conditions y{x, ) = f{x) and 
yt{x, 0) = g{x) is well known jl^. Ifsj] to correspond to d'Alembert's solutionlsi^: 



y{x, t) = - {fix + ct) + fix - ct)) + — / gis) ds (6) 



2 V / . V ''2c 



x+ct 



x~ct 



Hence initial displacements travel left and right. Initial velocity, however smears over a 
widening domain of influence. 
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Alternatively this results can be derived using the theory of partial differential opera- 
tors. Writing L = □ and X = [x, t) one arrives at the generic partial differential operator 
equation: 



LuiX) = k{X) (7) 

Here it is important to note that it is no longer required that u{-) is in the class of twice 
differentiable functions C^, but rather that u{-) and k{-) are distributions 32] or generalized 



'unctions. What this means in detail we will have to defer to expositions available elsewhere 
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loj l. For our purpose it is interest to note that in the theory of distributions jumps 
and discontinuity are gracefully and meaningfully included in the formalism including the 
definition derivatives of entities like the Dirac-delta 6{-) over the space of suitable functions 
testing for this property. In addition it makes continuous convolutions a central operation to 



the calculation of continuous solution, which make it's treatment very simi 



ar to the study 



of discrete models [l^l which is regularly used for digital waveguide models ^. 

In our treatment here we will closely follow Hulshof 'sllSl and Joshi and Wassermann's 

n n 

lecture notes jl8| and also Edwards' text[l^ which are more accessible than, for example the 
technical survey by Egorov, Komech and Shubin^^. The interested reader is pointed to 
the latter for statements of necessary theorems as well as proofs or detailed references to 
original proofs. 

The fundamental solution of the equation ((7j) the effect of the operator L on the distri- 
bution u when it sees as input a Dirac-delta (as noted earlier this is not a function, but a 

n 

distribution). In the digital signal processing literature [17| u is called the impulse response 
of L. In a dynamical sense it is the response to the inhomogeneous equation where the exter- 
nal force distribution is a Dirac-delta 6. Of immediate interest are solution forward in time 
the discussion is restricted to fundamental solutions in the positive half-plane with respect 
to time. This will be indicated by the superscript -|- to the symbol S for the fundamental 
solution: 

LS+{X) = 6{X) (8) 
It can be proved that the solution of both the homogeneous equation with initial value 



data and the inhomogeneous equation with some external force distribution can be recov- 
ered from the fundamental solution 
convolution in digital signal processing 
the force distribution and the initial value data u{x,0) = /(■) and Ut{x,0) = g{-). In the 



16 1 ■ T his is done, in analogy to the impulse response 
ig l31 by convolution of the fundamental solution with 



literature these are call Cauchy data. 

The fundamental solution 3, Q| of the one- dimensional wave equation can be derived 
to be: 



S+{x, t) = -H{t) [H{x + t)- H{x - t)] 



(9) 



Here H{-) is the Heaviside distribution, which is the distributional integral of the Dirac- 
delta distribution 6(-). In conventional functional form the Heaviside step- "function" can 

n 

be written as|19l|: 



H{x) 



X < 0, 

1 X > 0. 



(10) 



The interpretation of this equation is indeed important because it indicates, that the 
"system response" of a wave operator to an input impulse are not isolated traveling wave 
pulses but traveling step-distributions. 

It may be convenient to think of the fundamental solution as the "distributional continu- 
ous impulse response" . This makes sense because the solution of equation can be recovered 
by convolution of the fundamental solution with the Cauchy data. The continuous convolu- 
tion has the familiar form jisl . llG^ : 



u{x) = {S+ * f){x) 



S^{x — s)f{s) ds 



(11) 



where f{x) is Cauchy data in one variable. 

If both data and the fundamental solution are in two dimensions (as is possibly the case 
for external force distributions), one need to convolve over both variables. 
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Specifically it can be shown (see 13, Q) that for the set of Cauchy data u{x, 0) = f{x), 
Ut{x,0) = g{x) and Lu{x,t) = k{x,t) one gets the complete solution: 



u 



[x, t) = Uf{x, t) + Ug{x, t) + Uk{x^ t) (12) 



with the convolutions: 



Uf{x,t)=St{-^t)*f{x) (13) 
Ug{x,t)=8+{-,t)*g{x) (14) 
Uk{x, t) = £^ * k{x, t) (15) 

Performing the convolutions yields the the solution which is equivalent to the conventional 
inhomogeneous initial value solution of the wave equation jisl ]: 



y{x,t)=^ifix + t,0) + fix-t,0)) 

rx+t 

+ - g{s,0)ds (16) 



x~t 

t j'X+(t-T) 



H — / / k(s,T) ds dr 

2 Jo Jx-{t-T) 



For proofs of uniqueness see Egorov, Komech and Shubin ly] . In the absence of external 
forces this reduces to the initial value solution 

The theory of generalized functions for partial differential operators explains why equa- 
tions (0), derived for the forced case, and (jHI), derived for the homogeneous initial value case 
have similar structure. The inhomogeneous case can in a generalized sense be made to in- 
clude the homogeneous initial value problem (also called Cauchy problem). For example, if 
we symbolically write k{x,t) = 6{t)k{x,t) one sees that the external force response matches 
the response to the initial velocity. 



1 



x+t 



Ug{x,t) = - g{s,0)ds (17) 

•/ X — t 



1 



x+t 



u-^{x,t) = - / k{s,0)ds (18) 

^ Jx-t 
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Hence an external force distribution which is impulsive in time k{x, ■) is indistinguishable 
from the equivalent initial velocity distribution g{x). Conversely it is noteworthy that initial 
values do not prescribe a state of the string alone, but also prescribe a sudden onset of such 
state. Hence initial values are not necessarily a free-field solution of the wave equation, that 
is the state of a string in the absence of force. Rather impulsive onset states act like external 
forces. 

In particular one can write the solution of the wave equation for the inhomogeneous initial 
value problem in the following simple form jisl : 

u{x,t) = S+ * g{x) +S+ * f{x) + [ S+_^*k{x,s)ds (19) 

Jo 

Next we note the differentiation of distributions on step-functions (see also the Appendix 

El): 



(H', 0) = {5, 0) (20) 

observe, using this relationship, that the derivative of the fundamental solution Q are 
two propagating Dirac-deltas, here again in symbolic notation: 



g^+ = ^H{t){6{x + t) + 6{x-t)) (21) 

For the current discussion it remains only to point out that the first contribution of (fTU)) 
look like propagating impulses under differentiation. 

From (jl9|) we see that the solution of the wave equation will only stay on the characteristic 
lines ^ and t] for a restricted class solutions of the wave equation. Only for very exceptional 
cases of initial velocities g and external forces k will the solution not integrate into the 
domain. Rather generically one ought to expect them to integrate into the inside of the 
forward characteristic cone x ± t > as depicted in Figure ^ 

The condition in which integration into the interior of the characteristic cone does not 
occur will be discussed in the discrete case and in this context has been discovered by 
Karjalainen and Erkut 3|- For this to be appropriate, care needs to be taken when assuming 



the impulse response of the system to have a particular form, otherwise one ought to expect 



both contributions at the same time. Additionally, in a physically realistic situation, when 
it cannot be guaranteed that an excitation is of purely displacement-type, one ought to 
expect that the dynamics of the system to integrate and persist over the whole inside of 
the characteristic cone. It also is worthwhile to point out that this is the mathematically 
consistent solution of the wave equation [16^]. Hence, assuming that the wave equation is 
a reasonable model for a given physical situation, one ought to expect such a behavior to 
exist and be observable. 

As a note, we state, that it is in fact well known that the Huygens' Principle, the isolated 
propagation of wave fronts, only holds for d'Alembertians of odd spatial dimensions greater 



or equal to 3 
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, 1^ meaning that only in this case does the fundamental solution concentrate 
on the characteristic cone t'^ = Other cases, including one spatial dimension 13], have a 
wave influence inside the characteristic cone. The idea that the solution would concentrate 
on the characteristic cone for all dimensions was held for a long time by mathematicians 
working on the wave-equation until Hadamard opened the development and the situation is 
since well understood |2fl| . 

Regarding Huygens' principle, a short popular exposition can be found in whereas a 
comprehensive technical exposition can be found in |22| . 

III. COMPARISON OF LEAPFROG AND WAVEGUIDE SOLVERS 

Next the leapfrog finite difference scheme will be compared with the digital waveguide 
method j^. The full treatment of derivations will not be repeated here and the reader is 
referred to these sources for details. 

Now let us discuss the properties of the so-called leap-frog finite difference molecule for 
the wave equation. The explicit time-stepping equation reads 0, [l^ : 

y{n + 1, m) =y{n, m — 1) + y{n, m + 1) 

(22) 

— y{n — 1, m) 

where the relationship between the discrete time step T and the spatial discretization X is 
chosen to satisfy c = X/T. In this case the leapfrog molecule can be shown to be consistent at 
sampling points with the wave equation Il0| . It can also be shown that waveguide solutions 

Q n 

are solutions of the leapfrog |9|]. As Karjalainen observes, the converse does not hold p. 
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For future discussion we will use the following symbols: 



y+ = y{n + l,m) (23) 

y> = yin,m+l) (24) 

y< = yin.m- 1) (25) 

y_ =y{n-l,m) (26) 

The condition that an impulse at the root of the molecule will only create responses along 

the characteristics of the wave can be expressed by the condition ?/+ = 0, i.e. there is no 
data within the characteristic domain of the molecule. 

Hence the non-integrating molecule condition reads: 

y^+y^-y_ = (27) 
From this we get the relationship of waves on the characteristics to their sum: 

y<+y>= y- (28) 
The updating rules for waveguides js^ are: 

yi{n,m) =yi{n-l,m + l) (29) 

yrin, m) = yr{n — 1, m — 1) (30) 

with the external force rule: 



yi(n,m) = ^f{n,m) (31) 
yr{n,m) = l-f{n,m) (32) 



2 

The wave reconstruction rule is: 



y{n,m) = yi{n,m) + yr{n,m) (33) 
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in response to an external force function f{n,m). If we take an impulse of height at 
time n — 1 we get: 



yi{n, m-1) =yi{n-l,m) = f{n - 1, m) = -y_ (34) 

yr{n, m + 1) = yr{n-l,m) = f{n - 1, m) = ^?/_ (35) 

The reconstructed wave using (jHHj) is zero everywhere except at: 

y{n,m-l) =yi{n,m-l) = ]^y- (36) 

y{n,m+l) =yr{n,m + l) = ]^y- (37) 

y{n, m) = yi{n - 1, n) + yr{n - 1, m) = y_ (38) 



and we see that the non-integrating case of the leapfrog (pSj) is satisfied with y< = l/2?/_ 
and ?/> = l/2y_. 

The leapfrog will "integrate" whenever condition (PHj) is not satisfied. To study the 
behavior within the characteristic domain it is first assumed that the elements on the char- 
acteristic of the molecule ?/< and ?/> are unaltered. That is, the same waves as before travel 
outward in the molecule. This leaves us to study an altered relationship between ?/+ and ?/_. 

Let y- be the difference of , the molecule value for the non-integrating case (PH|l . and 
y_, an assumed contribution to the interior of the characteristic domain. Then we get: 

l/+=l/<+l/>-(y° +y_) (39) 
= y< + 2/> - 2/° (40) 

Subtracting (jiUI) from we get: 



y+ = y- (41) 

Hence the response at at the interior point of the characteristic domain is constant with 
regards to the contribution of the incoming wave that violates the non-integration condition 




To study how this behavior, one can illustrate the response of the leapfrog to an initial 

1: 



11111 
1111 

1 1 1 (42) 
1 1 
1 

and compare it to an excitation which observes (j2Hl)- y< = y> = ^ aiid = 2: 

1 1 
1 1 

1 1 (43) 

1 1 
2 

Observe that ()4H|1 appears to be the sum of traveling histories and they are time- 
symmetric around the intersection point. A time-symmetric solution is an equal contribution 
to the solution traveling forward and backward in time and their sum yielding the complete 
solution. 

D'Alembert's solution (jH)) can be used to investigate this observation when writing it in 
the form following Alpert, Greengard and Hagstrom 



yix, t) + yix, -t) =i {f{x + ct) + f{x - ct)) 

-j^ px+ct 

+ ^ / 9{s)ds 



x—ct 



+ ^(/(x-ct) + /(x + ct)) 



(44) 



nx—ct 

+ 9is)ds 

Jx+Ct 

Taking the time-symmetric sum we get: 



y{x, t) + y{x, -t) =f{x + ct) + f{x - ct) 



(45) 
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Similarly, by taking the difference, one finds the time-asymmetric case: 

px — ct 

y{x,t) - y{x,-t) = - gis)ds (46) 

C Jx+ct 

In the discrete case it is easy to see this property preserved in the leapfrog case: 

1 1 1 

1 1 
1 

(47) 
-1 

-1 -1 
-1 -1 -1 

1 1 

1 1 

2 (48) 
1 1 
1 1 

While (j^8|) nicely illustrates the time-symmetry and the "interference" of waves at the 
interaction point, ()47|1 is insightful as it clearly shows the properties of a velocity excitation. 
The displacement vanishes at the interaction moment, while the temporal slope is maximal. 
It should be made clear, that vanishing of data at one time-step in the leap-frog simulation 
does not imply velocity solutions. This can be seen if a positive and a negative impulsive 
wave cross, creating a time-symmetric situation that will not integrate inside the domain: 

-1 1 
-1 1 

(49) 
1 -1 
1 -1 

We observe that this situation does satisfy the time-symmetric equation (ji^ . 
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Alternatively similar results can be derived by discretizing the initial velocity g{x) directly 
using a matching center difference scheme [l^ : 



y+-y-=g (50) 
and including an arbitrary background field one gets: 

y^ = lir+n+9 (51) 

where we use the notation and to denote initial displacement wave contributions 
aligned with the left-right branch of the leapfrog-molecule. Observe that ()51|) does satisfy 
the same integration ()4H) and non-integration ()28|) conditions. Hence a velocity contribution 
g can be interpreted as any violation of the rule of the sum of incoming traveling waves. 

This behavior has been observed earlier. Karjalainen observed that an asymmetric pair 
of impulses need to be fed into a leapfrog motivated junction formulations to avoid in- 
tegration behavior 24]. The subsequent physical interpretation is derived in [71 from a 
center-difference time-discrete velocity excitation. An interpretation of this result follows 
next. 



IV. SINGULAR PROPAGATION FROM INTEGRATION 

The non-integrating condition can be algorithmically enforced by a method discovered 
by Karjalainen and Erkut 0,1^. Hence we will call this the Karj^alainen-Erkut condition. 
The rule is to present the excitation through a feed-forward filter 



a 



H{z) = l-z-^ (52) 
which can be derived from physical conditions by using a center difference velocity term 



We observe that the Karjalainen- Erkut condition ()52|1 creates two impulses from one and 
those impulses are center symmetric and sign-inverted. If we calculate those two impulse 
responses separately and then create the sum, we see that the impulsive propagating solution 
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comes about as the difference of two Heaviside distributions. The pulses are represented by 
their sign only as the amplitudes are assumed to be matched: 



+ + + + + 
+ + + + 

+ + + (53) 
+ + 
+ 

+ 

- - - 
- - 

0-0 (54) 





+ + 
+ 00 + 

+ + (55) 
+ + 
+ 

Hence we see that in an impulse-response interpretation of the leap-frog, a Heaviside 
integration over the characteristic cone of influence is sensible and the Karjalainen-Erkut 
condition ensures that each Heaviside integration is matched with a delayed sign-inverted 
response that cancels all interior integration of the first impulse to leave unaltered the 
traveling impulse solution. 

V. EFFECTS OF THE BOUNDARY 

Next the effect of imposing spatial boundary conditions is studied. For this it is assumed 
that the solution of the wave equation is only meaningful and defined for a compact domain 
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Q. The length of the domain is denoted by L = For the boundary of the domain we 
write dQ and the interior of the domain is defined by the quotient Q \ dQ. On each distinct 
point of the boundary dQ we impose one boundary condition. Fixed ends u{dQ) = we call 
Dirichlet boundary conditions, whereas open ends Ut{dQ) = we call Neumann boundary 
conditions. Note that a circular domain = 'u(O) is a periodic unbounded domain. 

The behavior at the boundary can be conveniently studied by extension of the domain. If 
the boundary is of Dirichlet type, the value of u needs to vanish at the boundary and hence 
the extension needs to be odd. In the case of Neumann conditions u needs to be even. As the 
resulting infinite extension is periodic in 2L this extension can be interpreted as a periodic 
unbounded domain of this length 13j. We will denote the original domain by subscript 
and extended domains by indicies n G Z \ 0. The periodicity implies Q^n + 2 = for all 
m G Z. 

The following discussion will be restricted to the behavior in response to the velocity term 
in (jni). Observe that with periodicity we can write the integral as the sum of contributions 
of the periodic domains: 



ri=x+t 



^=x-t 



g{s)ds = 



g{s) ds 



(56) 



Hence we integrate over all contributions above the characteristic lines from an excitation 
point of the periodic domain. 

With Dirichlet conditions one gets the odd extension 25|: 



g{x) = -g{2m\VL\ - x) (57) 
and for Neumann conditions we get the even extension: 



g{x) = g{2m\n\ - x) (58) 

with X G i^o S'lid m G Z \ 0. 

Integrating up to the point where the characteristic lines have reached 2\Vt\ one gets for 
Dirichlet boundary conditions: 
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u{Q) = / g{s) ds — 9{^^) ds 

Jno J^n±i 







(59) 



Hence integral contributions cancel every 2\^l\ and the maximum amplitude is bounded 
by the integral of g{-) over the original domain Qq. 

The same procedure for Neumann boundary conditions leads to: 



Figure El shows the behavior for a string tied at the ends (Dirichlet conditions) after an 
initial impulsive distribution. It reveals many properties of the effect of the boundary on 
the integration of velocities under Dirichlet boundary conditions. It shows the odd-periodic 
extension of the domain Qq to n G Z, it also shows the cancellation and constructive 
interference effect of overlapping integration regions. It also shows the 2\Q\ cancellation of 
waves. Erkut and Karjalainen[26| (compare their Figure 7) reported numerical simulations 
using the leapfrog molecule with comparable results, which hence matches the situation of 
the continuous model. 

A. Linear Growth of Displacement 

Observe that the Neumann condition leads to a linear increase in the displacement as a 
response to velocity or force data being present. 

The difference between the Dirichlet condition and the Neumann condition can be inter- 
preted as the difference between an alternating sum and an accumulative sum. 

In the Dirichlet case the sign of the area integrated over alternates with periodicity \Q\ 
and hence any finite bounded signal g{-) will produce an alternating sum which is bounded 
similarly but infinitely periodic. 

In the Neumann the signs match and hence the area integrated over increases with every 
iteration over the domain by the integral over the finite bounded signal g{-). Once the 
support of (?(■) has been exhausted, this obviously corresponds to a linear increase. 




(60) 
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This is however, not an unphysical situation. This corresponds to constant kinetic energy 
being present in the string and hence imphes that the energy is bound. This can easily be 
understood as hnear increasing displacement implies constant velocity, which in turn implies 
constant energy. Hence energy is conserved 3] • It can be interpreted as a string moving at 
constant velocity, which is meaningful as the Neumann conditions imply that the string is 
not tied down at the boundaries. Hence linear buildup in a displacement-like wave variable 

n 

is energy-conserving j35|]. 

Numerically this is still an undesirable situation because even in the absence of numerical 
imprecision, the dynamic range of numbers are bounded and hence an infinite increase cannot 
be represented. 

The case of Neumann boundary conditions is interesting because it highlights the differ- 
ence between notions of stability as customary in discrete signal literature [ij] and stability in 
physical situations. The Neumann displacement response to any bounded input will be un- 
bounded and hence is evidently not bounded-input bounded-output (BIBO) stable, see Op- 
penheim and Schafer[l3|, P- 20. We suggest that this BIBO-unstable but energy-conserving 
system be called physically stable. The B IB 0-inst ability is a discrete computational problem 
and not one 



VI. IMPLICATIONS 



This paper discussed the linear lossless wave equation and its relationship to discrete 
models, namely a finite difference scheme called leapfrog, and the digital waveguide method. 
It is shown that the waveguide model corresponds to the finite difference scheme in the 
absence of integration. In the continuous case, integration can be expected to occur when 
initial velocities or external forces are present. In this light the observed integrating behavior 
of finite difference discretization in the time domain using the leapfrog molecule displays 
results consistent with the continuous model. Here we assume that the wave equation is 
at least in principle physically meaningful for the modeled situation. If this is the case one 
should expect consistent behavior of the related discrete models. 

In relation to this argument, a use of waveguide discretization that does not include 
contributions inside the characteristic cone, does not include the integrating behavior of 
the model equation. In general both integrating and non-integrating responses are to be 
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expected and hence should be present unless they can be explicitly excluded for physical 
reasons. 

This also implies that the impulse response in just one variable in general will not carry the 
full dynamics of the string. Hence any assumption of the non-integrating impulse response 
in one variable in the construction of physical models might contain deviations for it only 
covers a reduced set of the solution space. 
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APPENDIX A: PROPERTIES OF DISTRIBUTIONS 

Let / be a distribution on a real open interval Q and let be in the the set of test 
functions l){fl) then one has [l^ : 

(0,/')= / f<P'df, = -{<!>' J) (Al) 
Jn 

and for arbitrary derivatives: 

(0,av) = (-1)1^1(5^0,/) (A2) 

the Dirac delta S has the property: 

(0,(5) = 0(0) (A3) 

hence returns the value of at 0. By the differentiation rule the higher order derivatives 
of the Dirac delta returns the higher order derivatives at with alternating sign: 
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= -(l)lPl9P0(O) (A4) 
Let H be the Heaviside distribution. It is defined as llSll: 



/CO POO 
H{s)^{s)ds= / (A5) 
-oo Jo 

It hence permits the positive part of (p over the domain. The derivative of the Heaviside 
distribution H yields (using (jAlj) and ()A4j) ) the Dirac-delta: 



(0,/f') = -(0',if) = 0(O) = (0,5) (A6) 

APPENDIX B: THE WAVE EQUATION AND FIRST ORDER SYSTEMS 

In order to derive the relationship between the wave equation to first order systems, we 
discuss two forms of such systems, namely, two transport equations in one variable and two 
transport equations in a mixed pair of variables. 

A generic version of a system of inhomogeneous first order hyperbolic equations reads: 



+ = h,{x) + hit) (Bl) 
ox at 

+ d% = h{x) + h,{t) (B2) 
ox ot 

For simplicity assume that the force terms are separated in the independent dimensions. 
Then a second order version is usually derived taking the derivative of one equation with 
respect to t and the other one with respect to x. The cross-term y^t can be eliminated and 
one gets two equations: 
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The key observation is that one second order equation ()B3|) or ()B4|) is not strictly equal 
to the system of first order equations ()B1|) and ()B2|) . It is only equal up to two functions 
(whichever got eliminated, hi,hi or /i2,^3)- They are equivalent up to two "constants of 
integration" . 

For systems of first order linear equations in two independent variables a related proof 
holds. An intuitive interpretation is that in fact for first order equations of the type: 



Ux + wt = giit) (B5) 
Wx + ut = g2ix) (B6) 

one sees that the reduction to second order equations in u by differentiating ()B5|) with 
respect to x and ()B6|) with respect to t one gets: 



Uxx -utt = (B7) 
-Wxx + wtt = gi(t) - g'^ix) (B8) 

Note that differentiation eliminated gi and g2 in one case and hence the homogeneous wave 
equation is again indistinguishable for both the homogeneous and a class of inhomogeneous 
systems of first order equations and in this sense they are equivalent only up to two functions. 
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Interior of the Characteristic Cone 




Characteristic Lines 



FIG. 1: The characteristic cone of the one-dimensional wave equation. 




FIG. 2: Leapfrog computational molecule for the one-dimensional wave equation. 




FIG. 3: Sum of velocity domains, ilo is the original string domain and r2„ with n G Z \ are 
domains created by continuation of the domain obeying the boundary condition u{dQ) = 0. 
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